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ABSTRACT 


The  purpose  of  this  research  is  to  study  the  Hough  transform  method,  as  applied  to 
the  detection  of  tracks  of  underwater  moving  objects  in  Lofargrams.  The  subjects  in¬ 
cluded  are  the  Hough  transformation,  clustering  study,  and  reconstruction.  Two  meth¬ 
ods,  LAS  cluster  technique  and  Sorting,  are  used  for  cluster  analysis.  Encouraging 
results  are  obtained  from  the  Sorting  method.  A  further  improvement  of  the  Sorting  is 
shown  to  yield  better  results  in  processing  noisy  track  data.  Experimental  results  dealing 
with  noise  free  artificial  data,  noisy  artificial  data,  and  real  noisy  data  are  presented. 
The  improved  Sorting  technique  presented  in  the  thesis  has  shown  improvements  com¬ 
pared  to  the  straight  forward  Sorting  when  it  is  applied  to  spectral  component  tracking. 
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I.  INTRODUCTION 


Track  detection  and  recognition  of  moving  objects  in  Lofargrams  is  a  subject  of  in¬ 
terest  in  image  processing.  Especially  in  military  application,  it  is  an  important  problem. 
The  spectral  tracks  in  a  lofargram  usually  exist  in  a  noisy  background.  To  detect  and 
recognize  a  curve  belonging  to  a  target  from  image  data  is  a  problem  of  importance. 
This  research  started  with  the  characterization  of  the  line  and  curve  representation,  and 
then  concentrated  on  studying  the  Hough  transform  algorithm. 

The  main  objective  of  this  study  is  to  search  for  a  method  that  can  be  used  by  a 
computer  rather  than  by  a  human  to  detect  the  tracks  of  an  object  and  to  extract  infor¬ 
mation  pertaining  to  the  image  data  source.  The  study  is  aimed  at  the  recognition  of  the 
tracks  of  underwater  moving  objects  with  possible  dynamic  Doppler  shifts. 

A  problem  of  interest  is  how  to  get  track  information  in  situations  where  the  object 
exhibits  dynamic.,  which  are  usually  observed  as  a  string  of  line  segments.  In  addition, 
how  to  determine  the  missing  data  (i.e.  temporary  disapperance  of  the  line)  from  the 
noisy  image  data  is  also  an  important  issue  to  be  studied.  These  will  be  discussed  as  a 
modification  of  the  Hough  transform  method. 

Before  experimenting  with  the  Hough  transform  method,  the  image  data  must  be  set 
up  in  the  right  format.  The  image  sources  can  be  real  data,  or  artificial  data.  In  the 
majority  of  experiments,  artificial  data  is  used.  For  further  study,  real  data  can  also  be 
used  in  the  computer  experiment. 

To  create  artificial  data,  the  user  can  either  write  his  own  program  or  use  the  com¬ 
mercial  signal  processing  software  package  called  Interactive  Laboratory  System  (ILS). 
The  tool  to  display  the  artificial  image  data  is  the  IM-4000  Image  Manager  in  the 
VAX, 'VMS  system.  The  IM-4000  Image  Manager  is  a  software  package  which  can  dis- 


play  images  up  to  an  array  size  of  512  x  512  pixels.  With  8-bit  pixel  resolution,  it  is 
possible  to  display  simultaneously  up  to  256  colors  or  gray  shades,  which  is  adequate  for 
the  display  of  most  digital  images.  A  functional  block  diagram  of  the  IM-4000  Image 
Manager  is  shown  in  Figure  1.  [Ref.  1] 

The  Hough  transform  method  allows  the  deteetion  of  line  structures  in  a  parameter 
domain.  This  method  is  stable  in  the  presence  of  noise.  Because  the  cluster  finding  de¬ 
termines  the  success  of  line  detection,  two  methods  were  used  in  the  study.  One  is  using 
the  routines  of  the  Land  Analysis  System  (LAS),  and  the  other  is  based  on  the  Sorting 
technique. 

In  order  to  use  the  routines  of  the  LAS,  a  splitting  technique  is  used  to  create  two 
bands  of  subimages.  One  is  for  p  band  and  the  other  is  for  0  band.  After  the  image 
splitting,  the  Land  Anabsi;  System  (LAS)  routines  are  used  to  find  the  position  of  the 
clusters.  In  this  study  HINDU,  KMEANS  ,  ISOCLASS  [Ref.  2]  routines  of  the  LAS 
are  adopted  to  do  the  cluster  analysis.  Comparisons  and  results  of  these  routines  are 
discussed  in  latter  chapters. 

To  confirm  the  accuracy  of  the  cluster  analysis  a  reconstruction  procedure  based  on 
the  positions  of  the  clusters  is  obtained  from  the  LAS.  If  the  Sorting  method  is  used,  a 
threshold  must  be  determined.  To  verify  the  performance  of  the  Hough  transform 
technique,  experiments  with  noisy  images  have  to  be  performed.  In  the  reconstruction 
with  noisy  images,  the  Hough  transform  method  can  find  straight  lines  easily.  For  dy¬ 
namic  lines  (i.e.  curves),  there  are  still  problems  in  the  recontruction  process  which  are 
due  to  the  simplicity  of  Sorting  technique.  Hence,  we  need  to  improve  the  Hough 
method  to  recognize  curves.  A  later  chapter  will  discuss  an  extended  Hough  method 
that  provides  satisfactory  results.  In  a  real  life  scenarios  the  information  obtained  from 
the  Hough  transform  method  can  be  used  in  a  strategic  decision  process  for  further 
analysis. 


Figure  1.  IM-4000  control  block  diagram  (adopted  from  [Ref.  I]). 
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There  are  five  chapters  in  the  thesis.  Chapter  I  is  a  introductory  description  of  the 
track  detection  and  the  line  representation  problem.  Chapter  II  presents  the  theory  be¬ 
hind  the  computer  simulation.  The  principle  of  the  Hough  transform  method  and  al¬ 
gorithms  for  curve  detection  will  be  discussed.  Chapter  III  emphasizes  the  detailed 
implementation  of  the  Hough  transform  method  such  as  clustering,  and  reconstruction 
in  noise.  Chapter  IV  shows  the  results  obtained  from  the  computer  experiments  in  the 
presence  of  noise.  An  improvement  to  obtain  satisfactory  results  is  discussed  as  well  as 
the  experiences  gained.  Chapter  V  presents  the  conclusions  and  recommendations  for 
future  study. 
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II.  LINE  AND  CURVE  DETECTION  USING  THE  HOUGH 


TRANSFORM  METHOD 

A.  TRACK  DETECTION  PROBLEM 

Acoustic  track  data  (i.e  images)  are  usually  displayed  as  intensity  variations  in  a 
hybrid  domain,  where  the  horizontal  axis  is  the  frequency  and  the  vertical  axis  is  the 
time.  The  received  acoustic  data  is  processed  to  extract  the  spectral  information  as  a 
function  of  time  (i.e.  waterfall  display)  and  presented  in  the  above  image  format.  A 
stationary  object  emitting  a  constant  sinusoidal  oscillation  will  appear  as  a  constant 
vertical  line  in  the  acoustic  image.  If  the  object  is  moving,  a  time  varying  Doppler  shift 
is  introduced  to  the  signal  and  the  tracks  will  become  curves  Under  this  scenario,  the 
track  detection  problem  becomes  a  special  line  detection  problem. 

B.  HOUGH  TRANSFORM 
1.  General  Statement 

In  the  field  of  digital  image  processing,  one  of  the  transform  techniques  that  can 
be  used  to  detect  lines  and  curves  is  called  the  Hough  transform.  From  a  mathematical 
viewpoint,  it  is  clear  that  the  shortest  distance  between  any  two  separated  points  is  a 
straight  line  with  the  equation  expressed  as  y=  mx+  c,  where  m  is  the  slope  and  c  is  the 
intercept  of  the  line.  In  this  respect,  the  line  can  be  assumed  to  be  a  set  of  points  with 
the  same  value  of  m  and  c.  The  Hough  transform  method  detects  lines  and  curves  based 
on  the  relationship  between  a  line  and  the  points  creating  it.  In  real  images,  a  blurred 
line  can  be  recognized  by  the  presence  of  a  group  of  colinear  or  almost  colinear  points 
in  a  parameter  space. 

In  general,  if  a  picture  contains  n  points  (i.e.  discrete  white  points  lying  on  a 
black  background),  then  there  are  —  possible  lines.  It  is  necessary  to  perform 
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- - ^  rt3  comparisons  to  determine  colinearity.  Hence  the  detection  of  straight 

lines  is  problem  in  computation  and  may  be  prohibitive  for  large  n.  In  view  of  the  above 
mentioned  difficulty,  a  method  to  change  the  original  colinearity  problem  to  one  that 
finds  concurrent  lines  was  proposed  by  Hough  [Ref.  3]. 

2.  Fundamental  Theory  and  Representation 

The  Hough  transform  is  a  mapping  between  the  coordinate  of  the  image  space 
and  the  coordinate  of  the  parameter  space.  The  transformation  equation  for  a  straight 
line  such  as  y  =  mx  +  c  in  the  image  space,  will  create  a  common  point  at  (m  ,  c)  in 
the  parameter  space.  A  problem  will  occur  when  the  line  in  the  image  space  is  a  vertical 
line  since  the  slope  is  undefined.  An  alternative  way  is  to  use  the  ( p ,  0)  parameter  space 
instead  of  (m  ,  e)  space.  In  the  ( p  ,  0 )  space  a  straight  line  becomes  xcos0  +  ysin0  = 
P- 

The  Hough  transform  can  be  used  to  deal  with  the  detection  of  specific  struc¬ 
tural  relationships  between  pixels  in  an  image.  In  a  picture  plane,  one  pixel  of  the  x-y 
space  can  correspond  to  a  sinusoidal  curve  in  an  (p  ,  0)  parameter  space  according  to, 

p  =  xcosQ  +  ysinO  (2.  i) 

Consider  the  image  space  in  Figure  2(a)  with  a  pixel  point  at  (x, ,  y,).  It  corre¬ 
sponds  to  a  curve  in  the  parameter  space  shown  in  Figure  2(b).  If  we  add  a  second  point 
(x2  >  yi)  into  the  same  image  space  as  shown  in  Figure  3(a),  we  can  get  another  corre¬ 
sponding  curve  in  Figure  3(b).  These  two  curves  intersect  at  a  common  point  of  (p  ,  0). 
The  fact  is  that  any  other  points  lying  on  the  line  formed  by  these  two  points  (x,  ,  y,) 
and  (x2  ,  y2)  will  correspond  to  a  set  of  curves,  and  these  curves  intersect  at  the  same 
common  point  (p' ,  0')  [Ref.  3].  Hence  for  line  detection,  as  long  as  we  can  figure  out 
the  location  of  the  common  cluster  point  in  the  parameter  space,  we  equivalently  de¬ 
tected  points  along  a  common  line  in  the  image  space. 
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Figure 


(a)  Image  space 


Image  space  and  parameter  space  transform 
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There  are  a  number  of  properties  for  the  Hough  transform  [Ref.  4].  These 
properties  are; 

•  A  point  in  the  image  space  (.v0)  y0)  correspond  to  a  sinusodial  curve  jrocos0  + 
>osin0  ~  p  in  the  parameter  space  (p  ,  9). 

•  A  point  in  the  parameter  space  (p  ,  9)  can  be  used  to  reconstruct  a  straight  line  in 
the  image  space. 

«  Points  lying  on  the  same  straight  line  in  the  image  space  correspond  to  sinusodial 
curves  touching  at  a  common  point  (p' ,  9 ')  in  the  parameter  space. 

•  Points  lying  on  the  same  sinusodial  curves  in  the  parameter  space  correspond  to 

lines  through  the  same  point  in  the  image  space. 

The  Hough  transform  can  be  used  successfully  to  detect  a  significant  distribution  in  the 
parameter  space.  In  the  presence  of  noise,  this  is  .omplicated  by  the  quantization  of 
both  the  image  (spatial)  space  axis  and  the  parameter  space  axis.  The  pixel  value  at  a 
spatial  position  (x  ,  y)  is  transfered  to  weight  the  sinusoidal  m  ipping  in  the  parameter 
space. 

C.  CLUSTER  ANALYSIS 

A  cluster  is  a  set  of  points  in  the  parameter  space  where  the  population  of  points  is 
high  compare  to  the  population  of  parameter  points  in  the  surrounding  region.  Locai 
ization  of  clusters  in  the  parameter  space  can  be  used  to  compress  image  data.  The 
quantization  procedure  is  an  important  step  m  processing  the  real  value  of  p.  The 
smaller  the  quantized  bin  size,  the  more  accuracy  can  be  obtained. 

In  this  research,  two  ways  are  implemented  to  find  the  positions  of  the  clusters.  One 
way  uses  the  software  routines  of  the  Land  Analysis  System  (LAS).  LAS  is  an  image 
analysis  system  designed  not  only  for  use  with  satellite  images  but  also  for  the  purpose 
of  manipulating  and  analyzing  digital  image  data.  It  includes  a  wide  range  of  functions 
and  statistical  tools.  The  second  way  uses  Sorting  for  cluster  analysis. 
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1.  LAS  Routines 

LAS  includes  a  variety  of  routines  for  cluster  analysis.  Three  routines  for  un¬ 
supervised  classification  (HINDU,  KMEANS,  and  ISOCLASS)  are  used.  They  are 
summarized  below. 

a.  The  HINDU  Classification  Routine 

HINDU  classifies  a  multiband  :m based  upon  its  multidimensional 
histogram  which  corresponds  to  the  disuHr'w^ri  Ir  J  *  parametric  space  in  our  study. 
Regions  in  the  histogram  with  high  density  art.  i  « mded  as  clusters.  The  us^r  specifies 
the  input  image,  the  minimum  and  maximum  acceptable  number  of  clusters,  and  the 
histogram  bin  size. 

b.  The  KMEANS  Classification  Routm? 

KMEANS  performs  an  unsupervise-..t  classification  using  the  K-means  al¬ 
gorithm.  The  basic  K-means  algorithm  <">’  crates  as  follows. 

Step  1:  Begin  with  an  arbitrary  set  of  cluster  center  for  the  desired  number  of  clusters. 

Step  2:  Compute  the  sample  mean  of  each  cluster. 

Step  3:  Reassign  each  sample  to  the  cluster  with  the  nearest  mean  distance. 

Step  4:  If  the  classification  of  ?-’l  samples  has  not  changed,  stop;  if  not,  go  to  step  2. 

c.  The  ISOCLASS  Classification  Routine 

ISOCLASS  performs  the  unsupervised  classification  of  an  image  using  an 
ISODATA-type  clustering  algorithm.  The  basic  ISCDATA  clustering  algorithm  oper¬ 
ates  as  follows 

Step  1:  Set  the  preset  count  for  printing  summary  cluster  statistics. 

Step  2:  Cluster  the  data  into  C  classes.  Eliminate  any  class  with  fewer  than  T  mem¬ 
bers. 

Step  3:  If  the  preset  count  has  been  reached,  go  to  step  6.  If  C  <  2N  split  any  cluster 

whose  samples  form  sufficiently  disjoint  g-oups.  If  any  clusters  have  her*,  opih,  go  to 

step  I. 
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Step  4:  Merge  any  dusters  whose  means  are  sufficiently  close. 
Step  5:  Go  to  step  i. 


Step  6:  Print  the  summary  cluster  statistics,  then  stop. 

In  the  algorithm  description,  C  is  the  number  of  classes,  T  is  the  minimum  number  of 
distributions  allowed  in  a  cluster,  and  N  is  the  approximate  desired  number  of  clusters. 

A  special  splitting  technique  was  adopted  when  we  use  LAS  to  tlnd  the 
clusters  in  the  parameter  space.  Because  uAS  accepts  band  images,  it  is  necessary  to 
split  the  c.iginal  image  into  sepera  ~  p  band  images  and  6  band  images.  Modification 
is  needed  in  processing  the  p  band  because  the  original  p  band  indue5;.  negative  values 
which  will  cause  errors  in  the  image  display. 

2.  Sorting  Technique 

Sorting  technique,  just  as  its  name  implies,  is  based  on  rearranging  the  distrib¬ 
ution.  Based  on  the  distribution  of  the  maximum  cluster,  the  user  can  set  up  a  threshold 
factor  to  pick  up  useful  clusters  that  are  sufficien.  large  enough  in  the  parametric  space. 
Threshold  is  ine\itabl>  used  in  both  the  LAS  and  the  Sorting  methods  to  decrease  the 
noise  effect  and  the  processing  time.  The  difference  is  in  the  sensitivity  of  the  threshold. 
Robust  threshold  is  usually  preferred.  Generally,  a  clustering  algorithm.is  based  on  it¬ 
erative  split  and  merge  operations.  msider  for  example  the  situation  that  each  cluster 
picked  up  from  Sorting  might  include  a  set  of  points  as  shown  in  Figure  4.  'f'hese  points 
can  result  in  a  worse  situation  in  reconstruction.  Hence,  it  is  necessary  to  use  a  simi¬ 
larity  measure  discrimination  to  merge  similar  points.  By  using  a  similarity  measure  in 
Figure  4  a  more  accurate  cluster  location  can  be  achieved  [Ref.  5:  pp.  419].  Figure  5 
shows  >\vo  duster  sites  v„  and  v,  given  by  v,=  (x„  y,)  and  iy  — (x;,jy),  respectively.  The 
simiiark  /  measure  between  v,  and  vy  is  defined  by 


<Vj,  Vj> 

5(V'  ’  V^~  <V[  ,  Vj>+<Vj  ,  Vj>-<Vt  ,  vy> 


(2.2) 


where 


<Vl,Vj>=\Vi\\fy  co s(v(,v)) 
=bxcxcosO 
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2 


and 


cosB— 


b2+c2-a2 
2  be 


b=j  (Xi-xof+fa-yo)2 

c=v  (xj-xo?+(jj-yo)2 

fl=\/  (xj-xf+tyj-yt)2 


and  (jr0,>-0)  is  the  origin  of  the  parameter  vector  space,  so 


S(vi ,  Vj)* 


<V,,Vi> 


<V(  ,  Vi>+<Vj  ,  Vj>- <V(  ,  Vj> 

bxcxcosd 


b2xc2— bxcxcosd 
b2+c2-a2 


,  *)  2  2 
b'+c  +a 


(2-3) 


The  similarity  measure  is  tested  ^ainst  a  threshold.  If  the  similarity  measure  is  less  than 
the  threshold,  v,  and  v}  sites  will  be  merged. 

D.  RECONSTRUCTION 

From  the  result  of  running  the  LAS  routines  or  the  Sorting  program,  the  positions 
of  the  clusters  are  known.  The  value  at  the  position  of  the  parameter  of  p  and  0  is  the 
accumulation  due  to  the  transform  and  the  track's  intensity.  Hence  the  reconstruction 
can  follow  a  reverse  procedure  as  follows. 
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For  a  given  image,  if  we  can  figure  out  the  cluster  in  the  parametric  space,  then  there 
is  no  problem  to  reconstruct  the  cluster  center  (p  ,  6)  by  following  the  equation 


A'  = 


p—y 


sind 

cosB 


(2.4) 


In  this  study,  we  will  select  a  reference  point  which  is  I  -rated  at  the  center  of  the 
image  space  by  shifting  an  amount  of.x0  and  j>0 ,  the  new  equation  is 

p  =  {x-  x 0)cos6  +  (y-  y0)sind  (2.5) 


and  the  new  reconstruction  equation  become 


*  =  *0  + 


p-{y-yo)sind 

cosd 


(2.6) 


E.  TRANSFORM  METHOD 

From  the  above  analysis,  a  computer  experiment  can  be  performed  step  by  step  in 
accordance  with  the  flow  diagram  in  Figure  6. 

The  procedure  consists  of  four  parts.  The  first  part  is  required  to  generate  an  image 
and  noise,  to  form  an  image  with  the  desired  signal  to  noise  ratio  (SNR).  A  program  for 
noisy  image  source  generation  is  listed  in  Appendix  A.  Also  an  image  with  a  given  SNR 
can  be  created  using  ILS.  The  operational  procedure  is  descriped  in  Appendix  B.  The 
second  part  is  the  Hough  transform  method  where  distribution  in  p,  6  parameter  space 
are  generated  and  quantized.  The  third  part  uses  one  of  two  different  methods  to  find 
the  positions  of  clusters.  The  fourth  part  is  to  reconstruct  for  experiment  verification 
the  image  from  the  selected  cluster  position  information. 

The  relation  between  image  space  and  parameter  space,  as  described  previously, 
suggests  the  following  algorithm  for  line  detection  [Ref.  6). 
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ARTIFICIAL 

IMAGE 


•  Table  look  up  technique  is  used  for  sinfi  and  cos0  calculation,  which  are  required 
in  the  computations  in  the  parameter  space. 

•  Select  the  center  of  the  image  space  as  the  reference  point. 

•  Compute  the  corresponding  p  value  by  using  equation  (2.5)  and  store  these  values 
in  a  two  dimensional  array. 

•  Quantize  the  p  values  and  6  values. 

•  Create  an  accumulator  array  jrt(0  ,  p) ,  whose  elements  are  initially  set  to  zero. 

•  Increment  all  quantized  points  in  the  accumulator  array  along  the  appropriate  line 
as  follows 

MO,  p)  =jri{6,  p)  +  ipd  .  (2.7) 

ic  =  ic+  1 

for  6  and  p  satisfying  equation(2.5).  Here  ic  is  used  as  a  counter,  ipd  is  the  gray 
level  of  the  pixel  that  is  added  to  each  element  in  the  accumulator  array. 

>•  Any  value  that  is  greater  than  the  threshold  located  in  the  accumulator  array  will 
correspond  to  the  most  likely  set  of  colinear  points  in  the  image  space. 

F.  NOISY  IMAGE  DATA 

In  the  ideal  case,  good  result  can  be  achieved  if  there  is  no  noise  involved.  But,  it 
is  inevitable  to  encounter  a  blurred  image  in  a  noisy  environment .  Hence  ,  how  to  de¬ 
crease  or  eliminate  the  effect  of  noise  to  recognize  the  track  is  an  important  problem  to 
be  studied. 

In  this  study,  artificial  image  data  is  generated  to  avoid  the  use  of  classified  infor¬ 
mation  as  well  as  to  have  truth  for  comparion  available.  The  noise  is  added  into  a 
noiseless  picture  to  provide  a  desired  SNR  [Ref.  7], 


17 


SNR(dB)  =  '  0  log 


(RMSS)2 

cx(RMSN)2 


(2.8) 


where  RMSS  is  the  root  mean  square  value  of  signal,  RMSN  is  the  root  mean  square 
value  of  the  noise.  For  a  given  SNR,  there  exists  a  proportional  constant  relationship 
between  SNR  and  c.  Hence  when  a  user  decides  to  create  a  SNR(dB)  image,  the  con¬ 
stant  c  can  be  computed  from  equation  (2.8).  Accounting  for  c  in  the  required  noise 
variance  to  obtain  the  desired  artificial  image,  is  straight  forward. 

The  principle  of  generating  a  desired  SNR  image  should  be  based  on  the  spectral 
domain.  That  is,  when  adding  the  noise  to  a  pure  signal,  it  makes  sense  to  add  signal 
in  the  spectral  domain,  which  is  equivalent  to  add  noise  to  the  signal  in  the  time  domain. 
The  additive  noise  is  Uniform  noise. 

The  generated  image  should  be  a  normalized  image  to  be  displayed  using  the 
IM-4000  software  package  for  verification  in  the  process.  Hence,  for  any  input  image 
it  is  necessary  to  test  the  gray  level  for  each  pixel.  The  value  of  the  gray  level  should 
be  restricted  to  an  integers  between  -128  and  127  as  a  signed  number  or  0  to  255  as  an 
unsigned  number  for  the  IM-4000  software. 


III.  IMPLEMENTATION  DETAILS  OF  THE  COMPUTER  SIMULATION 

In  this  chapter,  more  details  about  the  implementation  of  noise  and  signal  gener 
ation,  Hough  transform  method,  and  reconstruction  are  described. 

A.  SIGNAL  AND  NOISE  GENERATION 

According  to  ILS  tools  noise  is  created  randomly.  The  image  can  be  generated  in 
accordance  with  the  steps  in  appendix  B,  where  the  required  SNR  can  be  obtained  by 
adjusting  the  proportional  constant  c. 


RMSS 

nxln\0 

RMS*- 20" 


(3.1) 


B.  HOUGH  TRANSFORM  METHOD 

The  Hough  transform  method  is  applied  only  to  pixels  having  values  other  than  0 
in  the  image  array  which  requires  the  output  from  the  IM-4000  package  to  be  of  the 
unsigned  form.  The  reference  point  used  for  this  method  is  the  origin  of  the  coordinates, 
as  shown  in  Figure  7. 

Based  on  Figure  7,  the  range  of  p  becomes 


-rmax<Lp  <  rx 


max 


(3.2) 


where 


*max 


\fi* i 


2 

max 


i  -.2 

1  O'max 
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In  the  parameter  space,  the  p  axis  is  quantized  into  256  bins  while  Q  is  quantized 
into  180  bins.  However,  any  value  of  p  calculated  from  equation  (2.5)  is  a  real  number 
which  will  have  to  be  quantized  into  an  integer.  After  obtaining  the  integer,  the  negative 
p  values  still  need  to  be  mapped  into  positive  p  values  to  avoid  an  image  display  error. 
Besides,  the  positive  p  value  is  necessary  for  the  corresponding  quantized  accumulative 
array  to  have  an  accumulative  distribution.  In  order  to  meet  the  above  requirements, 
the  processing  is  done  in  the  following  fashion 

IR  =  (3.3) 

J2 

where  IR  not  only  stands  for  the  positive  p  value  but  also  stands  for  the  corresponding 
position  of  the  p  axis  in  the  parametric  space.  For  example,  take  an  artificial  image  in 
Figure  8(a),  after  executing  the  equation  (3.3)  every  p  value  transfered  from  the  image 
space  will  be  accounted  for  in  the  accumulative  array  to  form  an  image  array.  IM-4000 
can  display  this  image  array  as  shown  in  Figure  8(b).  The  more  often  a  given  position 
of  the  p  in  the  parameter  space  occurs,  the  bigger  the  value  of  the  accumulative  array 
and  the  brighter  the  gray  level  as  shown  in  Fig  8(b).  Any  value  of  the  element  of  the 
accumulative  array  greater  than  the  selected  threshold  will  represent  a  possible  cluster. 
For  illustration,  the  3-D  plot  of  Figure  9  shows  peaks  which  represent  the  possible 
clusters. 

The  next  step  is  to  find  out  the  positions  of  the  clusters.  Two  ways  are  used.  One 
is  to  use  the  LAS  routine,  and  the  other  is  to  use  Sorting. 

1.  Running  LAS  Routines 

HINDU,  KMEAN'S  and  ISOCLASS  are  three  different  routines  that  are  used 
to  find  the  sites  of  the  clusters.  Because  they  require  a  multiple  band  image  as  input,  it 
is  necessary  to  provide  a  normalized  p  band  image  and  a  6  band  image.  Two  2-D  arrays 
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(a)  An  artificial  image. 


) 


(b)  Result  after  Hough  transform. 
Figure  8.  Image  before  and  after  Hough  processing. 
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Figure  9.  3-D  plot  of  the  ( p  ,  0)  parametric  space. 


23 


are  used  to  store  the  value  of  p  as  well  as  the  the  value  of  6.  The  p  band  image  stores 
the  mapped  p  values  which  are  originally  obtained  from  equation  (2.5).  Before  the  p 
band  image  is  displayed,  the  negative  p  values  must  be  mapped  nto  positive  p  values  to 
avoid  an  error  both  in  future  cluster  localization  and  in  the  intake  band  display.  The  6 
band  image  stores  the  d  values  which  are  sequentially  assigned  between  0’  and  180’  for 
each  row.  For  the  example  of  Figure  8(a),  we  have  split  it  into  the  two  bands  as  shown 
in  Figure  10(a)  and  Figure  10(b).  While  running  the  LAS  routine  the  two  band  image 
data  will  contribute  to  the  cluster. 

2.  Normalization 

Normalization  is  a  necessary  step  to  transfer  the  negative  value  of  p  into  a 
positive  value  required  to  execute  LAS.  The  computation  follows  the  formula 

p'=(p  -  Prr.in)  X  factor]  (3.4) 


where 


factor  1  =  -3 — - 

rmax  —  rmin 

0'  =  (0  -  0min)  X  factor 2  (3.5) 


where 


factor!  =  -r - 3 - 

umax  -  ^min 

The  LAS  routine  will  accept  the  two  band  images,  then  use  merge  and  splitting 
techniques  to  eliminate  the  undesired  clusters.  The  result  is  statistical  information  about 
the  positions  of  the  clusters. 
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Figure  10.  Two  band  images  of  Figure  8. 
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3.  Sorting  Implementation 

There  is  a  direct  way  to  find  the  sites  of  the  maximum  magnitude  in  the  pa¬ 
rameter  space.  To  determine  the  clusters  a  threshold  is  necessary  in  this  approach.  In 
the  Sorting  procedure,  the  value  that  is  stored  in  the  first  address  of  the  jrt  array  will  be 
assumed  to  be  the  maximum  value.  It  is  compared  with  the  value  stored  in  the  next 
address  to  find  a  bigger  value  to  replace  the  previous  maximum  value.  Repeating  the 
comparison  throughout  the  whole  array,  the  peak  value  of  the  magnitude  in  the  pa- 
rametei  space  can  be  obtained  (i.e.  bubble  sort).  An  experimental  threshold  is  selected 
by  the  user  such  that  the  number  of  peaks  corresponds  to  the  number  of  lines.  Figure 
11  shows  the  Sorting  procedure  diagram. 

C.  RECONSTRUCTION 

The  value  of  the  positions  of  the  clusters  obtained  from  the  result  of  LAS  are  not 
the  final  information.  They  need  to  be  transfered  back  by  using  the  formula 


P  = 


foci  or  1 


+P; 


mm 


(3.6) 


■-£  .+<?  ■ 
factor 2  mm 


(3.7) 


If  we  apply  these  values  to  equation  (2.5)  and  equation  (2.6),  the  lines  can  be  re¬ 
constructed.  For  a  contenient  comparison,  the  example  in  Figure  8(a)  is  shown  in  Fig¬ 
ure  12(a)  again,  and  the  result  of  the  reconstruction  is  shown  in  Figure  12(b). 

In  this  chapter  miscellaneous  details  of  the  procedures  of  the  subroutines  have  been 
discussed.  It  provides  an  explanation  to  readers  wanting  to  execute  the  program.  The 
next  chapter  will  be  aimed  at  discussing  the  experimental  results. 
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(a)  Image  of  Figure  8(a). 


(b)  Reconstruction. 


Figure  12.  Image  example  reconstruction. 


IV.  EXPERIMENTAL  RESULT 


In  this  chapter  the  experimental  results  and  performances  of  different  approaches 
.re  presented.  These  experimental  results  reconfirm  the  potential  of  the  Hough  trans¬ 
form  method.  They  also  show  that  the  noise  has  less  effect  on  the  performance  when 
compared  to  other  image  detection  systems  such  as  the  Affine  Invariant  Matching  al¬ 
gorithm  [Ref.  8].  During  the  experimental  procedures,  the  cluster  analysis  is  done  using 
one  of  two  different  ways,  the  LAS  routines  or  Sorting.  In  using  the  LAS  routines,  there 
are  problems  that  can  not  be  solved  easily.  This  will  be  explained  at  the  end  of  this 
chapter.  Since  the  study  is  concentrated  on  the  track  analysis  of  the  moving  object,  the 
image  data  is  usually  a  segment  of  a  line  or  a  curve.  The  tracks  of  a  moving  object  w’ere 
generated  artificially  in  the  computer.  There  are  three  steps  in  the  execution  of  the  ex¬ 
periments.  First  we  work  with  the  noise  free  environments.  Then  we  work  with  the 
noisy  signal  to  test  the  tolerance  of  the  procedures  at  various  SXR's.  And  finally  we 
look  at  ways  to  improve  the  weakness  of  the  Sorting  method  to  obtain  more  reliable 
results. 

A.  EXPERIMENTAL  SORTING  RESULTS 
1.  Noise  Free  Tests 

First,  an  ideal  noise  free  situation  is  considered.  Artificial  images  with  three 
types  of  tracks  are  used.  These  types  are  lines,  curves,  and  mixed  curves. 
a.  Line  Test 

A  general  line  test  is  shown  in  chapter  II  as  an  example.  In  this  section,  a 
special  case  is  considered.  We  assume  that  the  unknown  artificial  image  has  two  vertical 
parallel  lines,  see  Figure  13(a).  According  to  the  Hough  transform  two  straight  lines 
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will  result  into  two  clusters.  But  the  results  show  four  clusters  as  seen  in  Figure  13(b). 
In  terms  of  the  3-D  plot,  the  corresponding  clusters  are  shown  in  Figure  13(c). 

The  reason  for  causing  four  clusters  is  that  the  cluster  points  will  be  located 
at  both  the  6  =  O'  and  the  d  =  180',  which  is  due  to  the  trigonometric  character  of  the 
transformation.  Two  cluster  points  exist  for  each  line.  One  is  at  6  =  0‘,  and  the  other 
occurs  at  6  =  179’.  In  reality  there  are  only  two  clusters  or  two  pairs  of  points.  During 
reconstruction  each  pair  should  reconstruct  its  original  line.  The  reconstruction  is 
shown  in  Figure  13(d)  proves  that  the  cluster  pair  explanation  is  correct. 
b.  Curves  Test 

In  this  section,  the  assumed  curve  image  has  a  shape  of  the  letter  S  as 
shown  in  Figure  14(a).  By  using  the  Hough  transform  the  clusters  are  not  easily  found 
from  Figure  14(b)  or  Figure  14(c).  The  result  of  the  reconstruction  is  shown  in  Figure 
14(d).  Though  the  Figure  14(d)  is  not  the  same  as  the  original  shape,  the  intersecting 
lines  more  or  less  provide  some  information  about  the  moving  track.  The  frequency  shift 
and  hence  differential  speed  of  a  moving  object  can  be  estimated  from  the  data  in  Figure 
14(d). 

c\  Mixed  Test 

A  S-curve  and  a  vertical  line  represent  the  mixed  test  and  are  shown  in 
Figure  15(a).  The  Hough  transform  yields  Figure  15(b)  and  the  3-D  plot  in  Figure  15(c). 
The  reconstructed  image  is  given  in  Figure  15(d).  The  reconstruction  is  not  identical  to 
the  original  so  that  some  modification  should  be  made. 

2.  Noisy  Test 

The  noisy  images  are  generated  by  using  the  ILS  software.  In  this  type  of  ex¬ 
periments,  signal  to  noise  ratios  are  varied  to  test  the  tolerance  of  the  procedures.  The 
tests  include  four  different  line  images  with  SNR's  of  20dB,  SdB,  3dB  and  -3dB. 
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(a)  Artificial  image(  straight  line). 


(b)  The  dense  region  indicate  the  clusters. 
Figure  13.  Line  detection  in  a  noise  free  situation. 
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(c!)  The  result  of  reconstruction. 
Figure  13  (continued). 
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(a)  Artificial  image(s-curve). 


(b)  The  dense  region  indicate  the  clusters. 
Figure  14.  Line  detection  in  a  noise  free  situation. 


(d)  The  result  of  reconstruction 


Figure  14  (continued). 
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(a)  Artificial  image  (s-curve  and  line). 


(b)  The  dense  region  indicate  the  clusters. 
Figure  15.  Mixed  line  detection  in  a  noise  free  situation. 
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(c)  3-D  plot  of  possible  clusters. 


(d)  The  result  of  reconstruction. 
Figure  15  (continued) 
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a.  Line  Test 


The  image  consists  of  five  vertical  lines  formed  in  ILS  with  five  frequencies  of  30,  50, 
60,  90,  120  Hz  respectively.  In  this  experiment  the  tracks  are  identified  in  the  scene  de¬ 
spite  the  noise.  Though,  a  filter  can  be  applied  to  remove  most  of  the  noise,  it  is  time 
consuming  in  running  filtering  while  processing  long  data  sequences.  The  Hough 
transform  can  solve  the  problems,  but  there  is  a  need  to  determine  the  threshold  to  pick 
up  a  sufficient  amount  of  clusters.  Figure  16(a)  shows  the  artificial  image  with  a  SNR 
of  20dB.  Using  the  factor  of  0.7  as  the  experimental  test  threshold,  the  test  recon¬ 
struction  recover  the  tracks  as  can  be  seen  in  Figure  16(b).  Note  the  line  corresponding 
to  the  lowest  frequency  has  been  split  in  the  reconstruction.  At  this  time  no  explanation 
for  this  phenomenon  has  been  found.  Also,  a  modification  of  the  algorithm  eliminates 
this  potential  problem.  Repeating  the  test  by  changing  the  SNR  to  8dB  is  shown  in 
Figure  17(a).  With  the  same  threshold,  the  reconstructed  results  are  shown  in  Figure 
17(b).  There  is  still  good  matching  between  the  signal  and  the  result  except  for  an 
artifact  of  a  swept  line.  The  test  image  for  3dB  and  -3dB  are  shown  in  Figure  18(a)  and 
lS(b).  The  reconstructed  images  are  shown  in  figure  18(c)  and  figure  18(d).  The  swept 
line  which  shows  up  in  Figure  18(c)  and  18(d)  is  due  to  noise  since  in  a  noise  free  simu¬ 
lation  no  swept  line  is  found.  These  results  can  be  improved  by  modifing  the  Hough 
transform  as  discussed  in  the  later  part  of  this  chapter. 
b.  Curve  Test 

The  S-curve  with  a  3  dB  SNR  is  shown  in  Figure  19(a).  The  reconstruction 
is  shown  in  Figure  19(b).  It  is  obvious  that  the  reconstructed  image  can  not  provide  the 
desired  information.  Hence,  for  curve  detection  it  is  necessary  to  make  a  modification 
in  the  Sorting  procedure  discussed  in  the  next  section. 
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(a)  Artificial  image  (line  with  SNR  20dB). 


(b)  Reconstruction. 

Figure  16.  Line  detection  with  SNR  of  20dB. 
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(a)  Artificial  image  (line  with  SNR  8dB). 


(b)  Reconstruction. 

Figure  17.  Line  detection  with  SNR  of  8dB. 
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(a)  Lines  with  SNR  of3d!L 


(b)  Lines  with  SNR  of-3dB. 
Figure  18.  Lines  detection  with  SNR  (a)  3dB  (b)  -3dB. 
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(c)  Reconstruction  of  3clB  image. 


(d)  Reconstruction  of  -3dB  image. 
Figure  18  (continued) 


(a)  Artificial  image  (curve  in  noise). 


(b)  Reconstruction  (original  approach). 
Figure  19.  Curve  detection  in  noise. 


42 


c.  Mixed  Test 

This  time,  a  Lofargram  is  taken  as  an  input  as  shown  in  Figure  20(a).  Here 
the  signal  to  noise  ratio  is  unknown.  The  Figure  20(b)  shows  the  result  of  the  Hough 
transform.  Again,  the  result  is  very  poor.  From  the  above  experiments  we  see  that  the 
noise  tolerance  of  the  line  detection  procedure  is  good,  while  the  noise  tolerance  for 
curve  detection  is  poor. 

B.  IMPROVEMENTS 

1.  Algorithm  Modification 

The  Hough  transform  method  has  difficulties  in  processing  curves  since  it  is 
difficult  to  pick  up  the  real  significant  clusters  without  picking  up  the  clusters  that  are 
caused  by  noise.  Hence,  the  first  thing  to  overcome  is  the  influence  of  noise.  One  simple 
step  is  to  split  the  whole  image  into  smaller  segments,  then  apply  the  Hough  transform 
to  each  segment.  Finally  we  put  each  iesult  to  its  corresponding  place  in  the  splitted 
image.  The  modified  algorithm  is  presented  as  following 

Step  1:  Split  the  whole  image  into  N  segment. 

Step  2:  Start  from  the  first  segment  and  compare  A  to  B,  where  A  is  the  mean  value 
of  the  segment,  B  is  the  mean  value  of  the  whole  image. 

Step  3:  If  A  >  B,  then,  run  the  Hough  transform  to  the  segment,  else  skip  to  analyze 
the  next  segment. 

Step  4:  Continue  step  3.  Stop  when  the  last  segment  had  been  compared. 

2.  Results 

The  Figure  21(b)  shows  the  improvment.  In  principle,  we  have  eliminated  the 
swept  line  caused  by  the  noise  in  Figure  21(a).  Figure  22(b)  shows  the  processed  result 
of  the  artificial  image  of  Figure  22(a). 

Figure  22(c)  shows  the  image  of  the  clusters.  Figure  22(d)  shows  the  reconstruction 
from  (c).  Figure  23(b)  shows  the  processed  result  of  the  artificial  image  from  Figure 
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(a)  Artificial  image  (mixed  in  noise). 


(b)  Reconstruction. 

Figure  20.  Mixed  curve  detection  in  noise. 
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(a)  Artificial  image  (same  as  Figure  18(b)). 
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(b)  Reconstruction  (after  modification). 
Figure  21.  Lines  defection  after  modification  from  Figure  18(b) 
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23(a).  Figure  23(c)  shows  the  Hough  transform  results.  Figure  23(d)  gives  the  recon¬ 
struction  from  (c). 

3.  Application  of  the  Result 

The  improvement  mentioned  above  shows  better  tracking  of  moving  object,  but 
multiple  lines  are  still  reconstructed.  What  needs  to  be  done,  is  to  analyze  the  recon¬ 
struction.  Take  the  Lofargram  shown  in  Figure  23(a)  as  an  example,  it  is  an  image  of 
the  time-frequency  space.  The  reconstruction  from  the  image  clusters  in  Figure  23(c) 
will  generate  the  corresponding  lines.  In  a  noisy  environment,  the  set  of  the  recon¬ 
structed  lines,  shown  in  Figure  23(d),  is  not  good  enough  for  further  analysis.  Here,  a 
repetition  of  the  Hough  transform  technique  is  suggested.  Take  the  image  in  Figure 
23(d)  as  an  input  image  to  a  second  Hough  transform  to  obtain  the  clusters  image  shown 
in  Figure  23(e).  By  using  a  clustering  method,  the  clusters  can  be  reduced  as  in  figure 
23(f).  Once  the  position  of  the  reduced  clusters  obtained,  we  can  implement  the  recon¬ 
struction  and  get  a  satisfactory  final  result  as  shown  in  Figure  24.  From  this  procedure 
slope  and  breakpoints  of  the  track  can  be  obtained  for  further  analysis. 

C.  EXPERIENCE  GAINED 

There  are  several  points  learned  from  the  experiments: 

1.  Unique  Coordinate  Reference 

The  coordinate  transformation  between  the  image  space  and  the  parameter 
space  is  important  not  only  before  cluster  analysis  for  accuracy,  but  also  after  cluster 
analysis  for  reconstrm  :ion.  Hence,  it  is  necessary  to  define  a  set  of  unique  coordinates 
to  avoid  errors.  For  example,  the  i  and  j  used  in  the  image  space  corresponds  to  6  and 
p  in  the  parameter  space,  and  will  be  used  throughout  the  program. 

2,  Boundary  Rule 

The  range  of  6  is  defined  as  0‘  ^  6  <  ISO'  by  Hough.  The  reason  to  check  ac¬ 
cording  to  this  definition  is  to  avoid  generating  cluster  pairs  in  an  ambiguious  fashion. 
This  was  discussed  in  connection  with  the  cluster  ambiguity  in  chapter  III. 
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(c)  The  image  of  the  clusters. 


(cl)  Reconstruction  from  (c). 
Figure  22  (continued) 
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(a)  Lofargram  image. 


(b)  The  reconstructed  Lofargram. 
Figure  23.  Curve  and  line  detection  on  a  Lofargram. 
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(c)  The  image  of  the  clusters  of  the  Lofargram. 


(d)  Reconstruction  from  (c). 
Figure  23  (continued) 
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(e)  The  clusters  from  (d)  by  a  second  Ilough  transform 
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(0  The  reduced  clusters  obtained  from  (e) 
Figure  23  (continued) 
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Figure  24.  Curve  and  line  detection  on  a  Lofargram:  (using  a  second  Hough 
transform). 
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3.  Normalization  and  Denormalization 

In  the  Hough  transform  negative  values  of  p  might  be  generated.  This  will 
cause  errors  in  finding  clusters  when  transfering  these  p  values  into  binary  format  to  be 
displayed.  Hence  the  p  negative  values  must  be  normalized  to  positive  values.  The 
computation  was  mentioned  in  chapter  II.  The  positions  of  the  clusters  obtained  from 
the  LAS  routine  provide  the  p  values  and  6  values,  respectively.  Since  they  will  not  be 
the  original  values,  inverse  processing  (denormalization)  in  accordance  with  the  equation 
(3.5)  and  equation  (3.6)  is  required  to  recover  the  original  p's  and  0's. 

4.  Quantization 

The  accumulation  array  is  quantized  along  the  p  and  6  axes.  For  instance,  assuming 
that  the  image  space  is  a  256  x  256  array,  then  the  accumulative  array  should  be  of  the 
same  dimension.  Basically,  the  larger  the  size  of  the  parameter  accumulation  array,  the 
more  accurate  the  procedure  will  be. 

5.  Experience  in  Using  LAS 

In  order  to  use  LAS,  the  splitting  technique  is  needed  to  split  the  input  image 
into  p  band  and  6  band  images.  But,  the  LAS  routines  have  difficulties  in  finding  the 
correct  positions  of  the  clusters. 

6.  Experience  in  Using  Sorting 

Sorting  technique  is  an  efficient  method  in  finding  clusters  by  setting  an  appro¬ 
priate  threshold.  But,  the  disadvantage  is  that  nearby  clusters  can  not  be  treated  as  one 
cluster.  This  is  the  inherent  problem  of  the  Sorting  technique.  To  solve  this  problem, 
cluster  analysis  needs  to  be  improved.  Of  course,  a  threshold  is  required  to  make  the 
similarity  decision.  From  experience,  using  a  similarity  measure  of  0.7  ~  0.9  yields  good 
result.  The  theoretical  consideration  for  an  optimal  similarity  measure  is  left  as  future 
work. 

7.  Noise  Tolerance 

The  Hough  transform  is  able  to  handle  noisy  signals.  The  signal  at  a  SNR  of 
-3dB  was  reconstructed. 
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V.  CONCLUSION  ■ 


The  purpose  of  this  research  was  to  extend  the  Hough  transform  method  to  the  de¬ 
tection  of  spectral  tracks  of  underwater  objects.  Data  used  were  artificial  images  gen¬ 
erated  by  computer  programs  and  Lofargrams.  Subjects  in  this  study  included  the  image 
space  transformation,  clustering  methods,  and  reconstruction. 

A  modified  Hough  transform  method  was  presented  to  yield  better  results  in  recon¬ 
struction  by  applying  segmentation  to  decrease  the  noise  effects.  The  clustering  was  an 
important  analysis  technique  to  search  for  significant  distributions  of  the  the  cumulative 
array.  This  step  significantly  influences  the  reconstruction.  Though  applying  the  LAS 
failed  to  find  the  accurate  cluster  sites,  an  alternative  uray  of  using  Sorting  can  make  it 
possible. 

The  track  detection  of  a  moving  object  is  an  interesting  problem,  and  is  of  interest 
in  military  applications.  This  thesis  discusses  a  practical  procedure  to  achieve  tracking 
analysis  of  moving  objects.  The  original  procedure  was  shown  not  to  work  well  with  real 
images.  An  improved  algorithm  yielded  better  results.  The  source  code  program  is  listed 
in  Appendix  C.  The  following  conclusions  can  be  drawn  from  this  research. 

•  The  Hough  transform  methodology  can  be  used  for  tracking  of  dynamic  lines.  This 
procedure  has  some  tolerance  to  noise. 

•  The  LAS  routines  did  not  provide  satisfactory'  results. 

•  The  Sorting  technique  provides  a  way  to  find  the  sites  of  the  clusters. 

•  The  Sorting  technique  has  difficulty  in  locating  a  single  cluster  point  per  track  in 
the  parameter  space,  and  hence  reconstruction  results  into  many  (false)  lines. 

•  A  threshold  factor  needs  to  be  determined  when  using  the  Sorting  technique. 
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APPENDIX  A.  NOISE  AND  SIGNAL  GENERATION  SOURCE  CODE 


This  appendix  contains  the  Fortran  source  codes  which  can  be  used  in  generating 
the  artificial  noise  and  signal  image.  Reader  can  create  the  desired  image  by  changing 


the  necessary  coefficients, 
program  sn-aen 

byte  a(256)/256*0/,y(256)/256*0/,b(256)/256*0/,\v(256)/256*0/ 
integer  m(l0),d(10),f(256,256),u(10) 
integer  h(256,256),max  1  ,min  1  ,ab,s,aa 
integer  ipd  1  ,summ, time, gray, sumf.dbvalue 
integer  mix(256,256),r(256),q(256),g(256,256), 
isy,s,aa,c,ic, ipd, sums 
real  ab  1  ,aa  1  ,f'ac.dbf,rmssl  ,adj(256,256), 
adjmix(256,256),y  1  (256) 
real  rmss.rmsn.db, index, ipd2 
open(unit=  l,name=  'sn256.dat',type=  'new', access 
1  =  direct'.recordsize  =  64) 
opcn(unit  =  2, name = 'n256.dat',  type = 'new', access = 

1  direct'.recordsize  =  64) 

opcn(unit=  3, name  =  's256.dat', type  =  'new', access  = 

1  direct'.recordsize  =  64) 

opcn(unit  =  4,  name = 'mix256.dat',  type  =  'new',  access  = 

1  direct'.recordsize  =64) 

byte  a(  1 28)/ 1 28*0/, v(  1 28)/ 1 2S*0/,b(I28)/ 1 28*0/, w(  128)/ 1 28*0/ 
integer  m(  10), d(  10), f(  128,1 28), u(  10) 
integer  ab,h(  1 28, 1 28), maxi  .mini 
real  abl,aal,fac 

real  dbfirmss  1  ,rmsn  1  ,adj(  128,1 28),adjmix( 1 28, 128), y I ( 1 2S) 

integer  ipd  1  ,summ, time, gray, sumf.dbvalue 

integer  mix(  1 2S,  1 28), r(  1 28), q(  1 28) 

integer  isy,g(!28,128),s,aa 

real  rmss,rmsn,db,index,ipd2 

integer  c,ic, ipd, sums 


£  +  A  V-  &  ?:■  <*•  »Ji  &  fy  fy  »*t  <:  i>  &  £  &  &  fe  %  ft  fa  j{s  jj:  »*:  >*t  £  »Jt  jjs  #  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft 

•  #!*.**  ^  .  _  »  *  _  _  f  A  A  y.  «  A  A  rf.  A  A  A  4.  A  A  A  A  4*  4*  A  A  A.  4.  V-  J.  A  4,  V.  J.  A.  V.  4.  4.  4.  A  V.  J.  4.  4.  A  rf.  4  •  A  A  A  A 

c  •  ••• crcat  signal  - *"■  ****•- . . . * ... * * ^ ... ... .. ... ... ... <. ... . * - ... ... ... ...  — - 

^  ft  ft 


dbvalue=  3 
m(l)=0 
m(2)  =  0 
d(!)  =  30 
d(2)  =  -20 
time  =  40 


gray=  100 
c=  2 
sums  =  0 
summ= 0 
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.  aa  =  0 

ic=0 
s=  11 
isy  =  128 
do  9  j=  l,c 
do  10  i=  l.isy 

fTi,j)  =  (m(j)*(i-64) + d{j))  +  64 

10  continue 

9  continue 

print  100,(U,l),fU/2),R2,l),H2,2) 

100  forma  t(  lx, 3(i5,lx),i5) 

do  11  j=  l,isy 
do  12  1=  l,c 
u(l)=llj,l) 

do  14  k=  l.isy 

if(k.eq.u(l))  g(k,j)=gray 

14  continue 

12  continue 

1 1  continue 

do  15  j=  l,isy 

do  16  i=  l,isy 
•pd  =  g(i,j) 
il(ipd.cq.gray)  then 
ic=  ic+  1 

sums®*  ipd**2  +  sums 
endif 

16  continue 

15  continue 
rmss=  sums**(0.5) 
print  17,ic,rmss 

17  format(lx,'ic  =  ',i5,3x,'rmss  =  ',f7.1) 
do  20  j  =  1  ,isy 

do  21  i=  l.isy 


a(i)=g(i.j) 

21 

continue 

20 

writc(l'j)  a 

continue 

c 

c 

c 


fc  #  >>  <5  <<  ❖  ❖  *  ❖  &  ❖  ❖  ❖  *  *  ❖  frffi  <«  #  <«  »>  t-  £  *  *  t-  <<  ❖  >C*  tr  £  ❖  fc  ❖  fc  *  ❖  <«  ❖  *  ❖  *  £ 

£  *  £  <*  crcftt  noise  ******,5t,^*<'^*******^****,c‘,>***5:,®‘ 

*  <*  #  #  *  »>  *  *  <S  £  v  *  *  *  *  *  *  »Js  *  ff  #  <s  #  »>  *3»  #  »>  *5t 


do  50  j=  l.isy 

do  51  i=  l.isy 

y(i)  =  (ran(s)-0.5)*tiine 

ab  =  y(i) 

aa  =  ab**2+aa 


h(i,j)  =  y(i) 

51  continue 

write(2'j)  y 
50  continue 


rmsn=aa*,>(0.5) 
print  52,rmsn 
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52  format(lx/rmsn  =  ',f7.2) 


c****  add  the  signal  with  the  noise  ********** 

£  %  jOuJs  i£«  lit  *  *  *  i>  &  *  *  <&  *  £  i>  £  *  «  Jjt  *  *  >2*  £  ft  »£  *  * £  Jfr  *  iji  *  *  *  #  *  *  *  *  *  i>  *  A  *  * 

do  60  j  =  l,isy 

do  61  i=  l,isy 

mix(ifj)=g(i,j)+h(i,j) 
ipdl  =mix(i,j) 

if(mix(i,j).gt.maxl)  maxi  =  mix(i,j) 
i^mix(i,j).lt.minl)  mini  =mix(i,j) 

61  continue 

60  continue 

fac = 255/(max  1  -min  1 ) 
do  62  j=  l,isy 
do  63  i=  l.isy 

mix(i,j)  =  (mix(i,j)-minl)*fac 
il{mix(i,j).gt,127)  then 
mix(i,j)=  mix(i,j)-255 
else 

mix(i,j)=mix(i,j) 

endif 

b(i)=  mix(i.j) 

63  continue 

writc(3'j)  b 

62  continue 
sm^summ^O.S) 
db= 20*LOG(sm/rmsn) 

index =  rmss/(rmsn*(EXP(.l  151293*dbvalue)-l)) 
print  53,db,index,sm 

53   rormatfJx/db  =  ',f7.2,3xf'indcx  =  ',f7.2,lx/sm=',f7.2) 


C** 


71 

70 

200 


redefine  db  to  assigned  DB 


sumf=0 
aal  =0. 
do  70  \  —  I,isy 

do  7!  i=  l.isy 

adj(i,j)= index*h(i  j) 
vl(i)=adj(i,j) 
abl  =yl(i) 
aal  =abl*lC‘2-}-aal 
continue 
continue 


rmsn^aal^O.S) 
print  200,  rmsnl 
format(  1  x/rmsnl  =  ',f7.2) 
do  73  j=  l,isy 

do  74  i=  l.isy 

adjmix(i,j)  =  g(i,j)+adj(i,j) 
ipd2=adjmix(i.j) 
sumf  =  ipd2**2 + sumf 
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if(adjmix(i,i).gt.l27)  then 
adjmix(i,j)  =  adjmix(i,j)-255 
else 

adjmix(i,j)=  adjmix(i,j) 
endif 

\v(i)  =  jnint(adjmix(i,j)) 

74  continue 
write(4'j)  w 

73  continue 

rmssl  =  suir<!-*(0.5) 
dbf=  20*LuG(rmssl/rmsnl) 
print  75, dbf, rmssl 

75  forma  t(  lx, 'dbf  =  ',f7.2,  lx, 'rmssl  =  ',17.2) 


APPENDIX  B.  COMMAND  OF  ILS  ARTIFICIAL  IMAGE  GENERATION 


In  this  appendix,  a  series  of  commands  will  be  presented  as  an  example  to  generate 
any  signal  with  appropriate  random  noise  as  follows. 

•  1:  SFIL  17 

•  2:  SCTX  128 

•  3:  STFU  SWD  1,256, 100, 256,,, 5 

•  4:  S1NA  SF256 

•  5:  SINA  II MY 

•  6:  SINA  N128 

•  7:  SSDI  LM  1.5 

•  8:  SFIL  17 

•  9:  SFIL  S19 

•  10:  SFIL  SI120 

•  11:  SNSI  Nl,256„8 

•  12:  SFIL  19 

•  13:  SFIL  S21 

•  14:  SCTX  256 

•  15:  SOPN  S 

•  16:  SSRE  1,128 
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17:  SDRE  SI, 3 


•  18:  SFIL  21 

•  19:  SFIL  S22 

•  20:  SOPN  S 

•  21:  SFFT  Pl,128„l 

•  22:  SFIL  22 

•  23:  DIO  RFN;WFRONT.dat 

•  24:  SR  SR2I 

From  1  to  7,  the  signal  will  be  generated.  From  8  to  11,  the  noise  will  be  generated. 
From  12  to  17,  the  mixed  signal  will  be  recorded.  From  18  to  24,  the  FFT  will  be  exe¬ 
cuted.  The  last  step  is  to  transfer  ASCII  code  into  binary  code  and  normalize  for  the 
display. 
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APPENDIX  C.  MAIN  SOURCE  CODE  PROGRAM 


program  mod5 
byte  a(256),p(256) 
integer  zz(16),split/64/,splitl/4/ 
real  threshold/. 95/ 

integer  k,m,c,u(  1 6),max,max  1  ,ta(256),tc(256),x  1(16,16) 
integer  max2,min2,max3,isy  1/256/ 
integer  ip(256,256),jrt(256,256),io,it,ir,ipd,td(256) 
real  mean.meanl 

integer  adx,ady,adt,adu,ip  1(16,1 6),jrt  1  ( 1 6, 1 6) 
real  tb(256),x0,y0,row(256),x(  16,16) 
real  sc(  1 6,2),dth,th,rowm,ro\v0,drovv,factor 
integer*4  isx/4/,isy/4/,ith/4/,iro/4/ 

open(unit  =  1  .name  =  'bcam00t.dat',  status  =  'old', access  =  'direct' 

1  ,recordsize=  64) 

open(unit=  3,name=  'tl283.dat',  status"  'new',  access  =  'direct' 

1  ,recordsize=64) 

open(unit  =  1 1  .name  =  'rcbeam00t.dat', status  =  'new', access  =  'direct' 

1  ,recordsize  =  64) 

c  ***»read  in  the  image  data  &  convert  to  integer.*******-********** 
sum  =  0 
icc=  0 

do  10  j=  l.isyl 
rcad(l'j)  a 
do  20  i=  l.isyl 
ip(i,j)  =  a(i) 
jrt(i,j)  —  0 

20  continue 

10  continue 

do  21  j=  l.isyl 
do  22  i=  l.isyl 

if(ip(i,j).lt.O)  ip(i,j)  =  ip(i,j)  +  256 

22  continue 

21  continue 

do  23  j  =  l.isyl 
do  24  i=  l.isyl 

sum  =  sum  +  ip(i,j) 

24  continue 

23  continue 

mean=  sum/(isy  l^isyl) +  20 
c  print  801, mean 

c801  format(lx,'mcan  =  ',f5.1) 
do  25  jj=  1, split 

do  25  ii=  1, split 
suml  =  0 
do  26  1=  l.splitl 

do  26  k=  l.splitl 

suml  =  suml  +  ip((ii-l)*splitl  +  k,(jj-I)*split  1  + 1) 
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,  ipl(k,l)=  ip((ii-l)*splitl  +  k,(jj-l)*splitl  +  1) 

26  continue 

icc  =  icc+  1 

mean  I  =  suml/(splitr**2) 
c  print  800, icc, mean  1 
c800  format(lx,'icc=  ',i6,2x,'meanl  =  ',f6.1) 
if(meanl.ge.mean)  then 

c  *  *  *  *  *  *  *  *  *  *  *  ^  0  u  g  j  ^  t  ra  n  s  fo  r  m  *  *  *  *  *  *  *  *  *  *  *  *  *  *  ***********  * 

call  ht2(ip)jrt,sc,iro,isx)isy,ith,drow,row0,x0)y0,ipl,jrtl) 

„*****  *  *  *  *  *  *  *  *  *  *  *  *  •**  * *  *  *  *  *  *  *  «  *  *  *  *  *  *  *  *  *  *  * .is  **********  *  *  *  * 


max  =  jrtl(  1,1) 
min= jrtl(l,l) 
do  30  j  =  1  ,isy 
do40i=l,isy 

ifl[]rtl(i,j)  .gt.max)  max= jrtl(i.j) 
if(jrtl(i,j)  .It.  min)  min  =  jrtl(ij) 
ta(i)  =  jrtl(i,j) 

40  continue 
30  continue 

pai  =  3.1415926 
c=0 

maxi  =  jnint(threshold*max) 
do  90  j  =  l,isy 

do  100  i=  l,isy 

if(jrtl(i,j)  .gt.  maxi)  then 
c  =  c  -l- 1 

tb(c)  =  (i-  r)*l  80/isy 
tc(c)=j 
endif 

100  continue 

90  continue 
do  1 10  j=  l,c 

ro\v(j)  =  tc(j)*drow-rowO 
110  continue 
do  120  j=  l,c 
do  130  i=  l,isy 

x(i,j)  =  x0+  (row(j)-(y0-i)’5'sin(tb(j)’'‘'pai/l  80)/cos(tb(j)’:'pai/l  80)) 
xl(i,j)  =  jnint(x(i,j)) 
td(i)  =  xl(i,j) 

130  continue 
120  continue 

do  140  j=  l,isy 

do  150  1=1, c 
u(l)  =  xl(j,l) 

150  continue 

do  160  m=  l,c 

do  161  k=  l,isy 
if(k.eq.u(m))  xl(k,j)=255 
161  continue 

160  continue 

140  continue  , 
do  162  j=  l,isy 
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do  163  i=  l,isy 

il(xl(i,j).ne.255)  then 
xl(i,j)=0 
else 

xl(i,j)  =  xl(ij) 

endif 

jrt((ii-l)*splitl  +  i,(jj-l)*splitl  +j)  =  xl(i,j) 

zz(i)  —  jrt((ii-l)*splitl  +  i,(jj-l)*splitl  +j) 

163  continue 

162  continue 

endif 

25  continue 

do  70  j=  l,isyl 

do  80  i=  l.isyl 

if(jrt(i,j)  .ge.  127)  then 
jrt(i,j)  =  jrt(i,j)-256 
else 

jrt(i,j)= jrt(i,j) 
endif 

a(i)= jrt(i.j) 

80  continue 

write(H'j)  a 

70  continue 

Q  *  *  J{!  *  *  *  *  *  *  *  *  *  *  *  *  *  *  .1!  *  *  *  *  *  *  *  ******  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  * 

call  modi 

£******  lit  «!•  v  ij!  ijj  v  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  ******  *  *  *  *  *  *  * 

call  t2563 

Q  *!•  &  v  v  -l'  -1"  v  v  v  ^  -I-  »**  v  v  t  v  >■{!  v  <{!  v  -Inj!  v  ^  if  I*,!  v  >1‘  v  ■!’  <1’  ^  >!> 

call  xx  • 

q**:  '••v  •»»  v  v  -r  «i*  v  v  •  -*•  ■('  v  >i>  v  v  4-  >!"  *  >i*  <i>  v  'i"  v  -I-  v  v  Jl*'!’'  i  >!'  •I"  -I-  >{•  v  *  •I’  'i*'  •1'  •!'  <!’•  4*  •<-  *  "i-  -I1  v  4  v  *  il*  v  *r  ■-i>  *i*  •!• 

call  clustering 

Q  *!•  »#*  v  V  *i*  V  V  V  V  V  V  V  V  2*.  V  *|t  V  ❖  *  *  *  *  £  *  ❖  *  <•  *  *  *  *  *  .1.  tf  *  v  v  v  *r  >1’  v  'i'*  ’r  if  -i-  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  * 

end 

.  V.  V.  J.  A  4  J  4  4.  4  4.  4.  4.  4.  4  4  4.  4.  4.  4  4.  4.  4*  4±  4.  4*  4*  4.  4.  4-  4.  4.  4.  4-  4.  4.  4.  4  4.  4.  4  f.  4-  4.  4,  4-  4,  4.  4-  4 .  4.  4-  4.  4.  4.  4 4.  4.  4.  4.  4 .  4.  4.  4.  4.  4 .  4-  4.  4.  4. 

A  >r  . . .  »»•  •  -t'  •f  't*  •>  v  *>•  >c  >i>  'i<  -r  *•*  ■>  »»■  •!••>>  v>#*v  >1*  *»•»■«•  »i*  >r  »*»•*•  ■>  'i*  »*-  9vv  3l*  '(•  •>  »»•  >r  -i"  •••  •>  >r 

q  >f  if  if  if  >1'.  i,v.  jjs  j);  jJ:  ;.v.  >.H  )).  if  »Js  >}t  >Js  if  >!<  »Ji  >Jt  i).  &  >Js  £  i).  .*}t  )J-.  if  if  if  if.  if  if  if  if.  if.  if.  if  if  >y.  )f  if  jJ;  <<  >Js  if  if  if  /f  if  if  if  if  if  if  if  if 


subroutine  modi 

byte  a(256),p(256) 

integer  zz(  1 6), split/ 1 28/, split  1  /2/ 

real  threshold/.95/ 

integer  k,m,c,u(  1 6), max, maxi ,ta(256),tc(256),xl ( 1 6, 1 6) 
integer  max2,min2,max3,isy  1/256/ 
integer  ip(256,256),jrt(256,256),io,it,ir,ipd,td(256) 
real  mean, mean  1 

integer  adx,ady,adt,adu,ip  1  ( 1 6, 1 6),jrt  1(16,16) 
real  tb(256),x0,y0,row(256),x(  16,16) 
real  sc(16,2),dth,th,rowm,row0,drow,factor 
integer’^  isx/2/,isy/2/,ith/2/,iro/2/ 

open(unit=  11, name  =  'rcbeamOOt. dat', status  =  'old',  access  =  'direct' 

1  ,rccordsize  =  64) 

open(unit  =  1 2, name  =  'rcbml  .dat',  status  =  'new', access  =  'direct' 

1  ,rccordsize=  64) 

c  >»*<'*read  in  the  image  data  &  convert  to  integer.********®******** 


64 


r' 


20 

10 


22 

21 


24 

23 


sum=0 
icc=0 

do  10  j=  l.isyl 
read(l'j)  a 
do  20  i=  l.isyl 
ip(i.i)  “  a(0 

jrt(i.j) = 0 
continue 
continue 
do  21  j=  l.isyl 

do  22  i=  l.isyl  , 

if(ip(i,j).lt.O)  ip(i,l)  =  ip0.l)  +  256 
continue 
continue 
do  23  j=  l.isyl 
do  24  i=  l.isyl 

sum=  sum+ip(i,j) 
continue 
continue 

mean  -  sum/(isyl*isy  l)+25 
do  25  jj  =  1, split 

do  25  ii=  1, split 
suml  =  0 
do  26  l—  1  .splitl 

suml  =  suml  +  i  p  ( ( i  i  - 1 )  ^  s  pH  j  j  -1)  ^  spl  1 1 1  + 1) 

ipl(.k,l)  =  ip((ii-l)*sphtl  -P  k,<jj- 1)  •  splitl-P  1) 
continue 

icc  icc4  1  meanlBSuml/(splitl**2) 
ifTmcanl  .ge.mean)  ft® 

c***«.*******hoUgh  transform- •  n 

call  Iu2(ip,irt,sc,iro,isx,isyjth,djro\v,ro\^,x0j^,jplt,j^j!()!i!# 

J,.  £  Jj.  J>.  .J.  A  J;  jJ*  f-  ^  *  *  /.{  $  *  *  $  *  **•  •"  *•'  . 

max= jrtl(l.l) 
min®  jrtl(l.l) 
do  30  j=  l.isy 
do  40  i  =  1  ,isy  . 

illjrtl(i.j)  .gt.max)  max  =  jrtl(i.j) 
if(jrtl(i,j)  -It.  min)  min  =  jrtl(i.j) 
ta(i)= jrtl(i.j) 
continue 
continue 
pai  =  3.141 5926 
c  =  0 

maxi  =  jnint(threshold*max) 
do  90  j=  l.isv 

do  100  i=  l.isy 

ifCjrtl(i.j)  -gt.  maxi)  then 
c  =  c+  1 

tb(c)  =  (i- 1)*  1 80/isy 
tc(c)=j 


26 


40 

30 
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endif 

100  continue 

90  continue 
do  110  j=  l,c 

row(j)=  tc(j)*drow-rowO 
110  continue 
do  120  j=  1  ,c 

do  130  i=  l,isy 

x(i,j)  =  xO  +  (row(j)-(yO-i)!i!sin(  tb(  jj^pai/ 1 80)/cos(tb(j)*pai/ 1 80)) 
x  1  (i,j) = jnint(x(i,j)) 
td(i)  =  xl(i,j) 

130  continue 
120  continue 

do  140  j  =  l,isy 
do  150  1=  l.c 
u(l)  =  xl(j,l) 

150  continue 

do  160  m=  l,c 

do  161  k=  l.isy 
if(k.cq.u(m))  xl(k,j)  =  255 
161  continue 

160  continue 

140  continue 

do  162  j=  l,isy 
do  163  i=  l,isy 

if(xl(i,j).ne.255)  then 
xl(i,j)  =  0 
else 

xl(i,j)  =  xl(i,j) 

endif 

jrt((ii- 1  )*split  I  +  i,(jj- 1  )*splil  1  +  j)  =  x  l(i,j) 

zz(i)=jrt((ii-l)*splitl  +  i,(jj-l)*splitl  +j) 


163 

continue 

162 

continue 

endif 

25 

continue 

do  70  j  =  l.isyl 

do  80  i=  l.isyl 

if(jrt(i,j)  .ge.  127)  then 
jrt(i,j)  =  jrt(i,j)-256 
else 

jrt(i,j)  =  jrt(i,j) 
endif 

a(i)  =  jrt(i,j) 

80  continue 

\vrite(I2'j)  a 

70  continue 
end 

q#  #  £  £  &  &  #  £  &  ft  >!<  Jqj.  JoilblC  call  *££***>**fc<5****,Ct,!,*<,5!#,5!*,C,*****,i!**,®t*,Ss**!*,5',fi‘ 

£  if.  *  <» £  « &  <i  <*  »>  *  £  *  f?  *  <*  *  £  *  *  #  to  *  #  ❖  ❖  *  <*  #  ❖  #  >'/  if-  *  *  <<  *  *  *  *  jfle  ijt *  # ijc  *  *  *  £  *  *  #  *  #  #  <*  *  ❖  *  £  #  £  *  <*  ❖  ❖ 


66 


subroutine  ht2(ip,jrt,sc,isx)isy,ith,iro,drow,row0,x0,y0,ipl,jrtl) 

y-  *  ■»,  v.  y.  y>  y.  y.  y.  y.  a.  y.  y  y.  y.  y.  y.  y.  y.  y.  y.  y.  y.  y.  y.y  y.  y.  y.  y  y.  y.  y.  a  y.  y.  y.  y.  y_  y.  y.  y.  y.  y.  y-  y-  y.  y.  y.  y.  -i.  y.  y»  y.y>y.y.  y,  y.  y.y.y.y,y.y, 


dimension  ipl(isx,isy) 
dimension  jrtl(ith,iro) 
real  sc(ith,2) 
real  xO,yO 
real  pai 

pai=  3.1415926 
jsx  =  isx 
jsy  =  isy 
jth=ith 
jro  =  iro 
do  100  j=  l,jro 
do  100  i=  l.jth 
jrtl(i,j)  =  0 
100  continue 

dill  =  pai/float(jth) 
do  10  i=  l.jth 
th  =  float(i-l)*dth 
sc(i,l)  =  cos(th) 
sc(i,2)  =  sin(th) 

10  continue 

rowm=  sqrt(fioat(jsx*jsx)  +  float(jsy*jsy)) 
row0=  rowiii/2.0 
drow=  rowm/fioat(jro) 
x0=  Iloat((jsx  +  2)/2) 

>0  =■ •flo<u((jsy  +  2)12) 
do  40  iy=  l,jsy 
>'=>>' 

do  30  ix=  l.jsx 
ipd  —  ipl(ix.iy) 
ii'(ipd.lc.0)  go  to  30 
x  =  ix 

do  20  it  =  l.jth 

row = (x-x0)*sc(it.  1)  +  (y0-y)*sc(it,2) 
r=(ro\v+  ro\vO)/drow 
ir = jnint(r) 

jrt  1  (it.ir)  =  jrt  l(it.ir)  +  ipd 
20  continue 

30  continue 
40  continue 
return 
end 

>>£$&&&£<<$&  4  $  &  &  &  &  £  &  &  *  £  &  $  ft  44  ft  £  &  &  £  tt. 


subroutine  t2563 
byte  a(256),p(256) 

integer  k,m,c,u(256),max,maxl  ,ta(256),tc(256),xl(256, 256) 
integer  max2.min2.max3 

integer  ip(256,256),jrt(256,256),io,it,ir,ipd,td(256) 

real  tb(256),x0,y0,row(256),x(256,256) 

real  sc(25612),dth,th,rowm,rou;0)drow, factor, threshold 
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integer*4  isx/256/,isy/256/,ith/256/,iro/256/ 

open(unit=  12,name=  'rcbm  l.dat', status  =  'old',  access  =  'direct' 

1  ,recordsize  =  64) 

open(unit=  1 3, name  =  'rcbmlast.dat', status  =  'iiew',access  =  'direct' 

1  ,recordsize  =  64) 

c  ****read  in  the  image  data  &  convert  to  integer.***********'***** 
threshold  =  0.32 
dol0j=l,isy 
read(12'j)  a 
do  20  i=  1  ,isy 
ip(i,j)  =  a(i) 

20  continue 

10  continue 

do  21  j=  l,isy 
do  22  i=  l,isy 

'  if(ip(i,j).lt.O)  ip(i,j)=ip(i,j)+256 
22  continue 

21  continue 

^  *  *  *  *  transform************************** 

call  htl(ip,jrt, sc, iro.isx.isy.ith.drow.rowO.xO.yO) 

niax= jrt(  1,1) 
min  =  jrt(  1,1) 
do  30  j=  l.isy 
do  40  i  =  1  ,isy 

if(jrt(i,j)  -gt.max)  max=jrt(i,j) 
il(jrt(i,j)  .It.  min)  min  =  jrt(i,j) 
ta(i)  =  irt(i,j) 

40  continue 
30  continue 

pai  =  3. 1415926 
c=0 

maxi  =  jnint(threshold*max) 
print  8 1  ,max 

81  format(lx,'max.  =  ',i5) 
print  82, max  1 

82  format(lx,'maxl  =  ',i5) 
do  90  j=  l,isy 

do  100  i=  l,isy 

il(jrt(i,j)  .gt.  maxi)  then 

c  =  c+  I 


tb(c)  =  (i-l)*180/isy 

tc(c)  =  j 

endif 

100  continue 

90  continue 
print  91,c 

91  format(  lx, 'counter  =  ',i5) 
do  110  j=  l,c 

row(j)  =  tc(j)*drow-row0 

110  continue 

print  92,ro\v(l) 
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>92  format(lx/rowl  =  ',f6.1) 

print  93,row(2) 

93  format(lx/row2==,,f6.1) 

do  120  j=  l,c 
do  130  i=  l,isy 

x(i,j)  =  xO  +  (row(j)-(yO-i)*sin(tb(j)*pai/l  80)/cos(tb(j)*pai/l  80)) 
xl(i,j)  =  jnint(x(i,j)) 
td(i)  =  xl(i,j) 

130  continue 
120  continue 

do  140  j=  l,isy 
do  150  1=  l,c 
u(l)=xl(j,l) 

150  continue 

do  1 60  m  =  1  ,c 
do  161  k=  l,isy 

il(k.cq.u(m))  xl(k,j)  =  255 

161  continue 

160  continue 

140  continue 

print  151,u(l),u(2) 

151  rormat(lx,i5,lx,i5) 
do  162  j=  l,isy 

do  163  i=  l,isy 

if(xl(i,j).nc.255)  xl(i,j)  =  0 
163  continue 

162  continue 

do  70  j  =  1  ,isy 
do  80  i=  l,isy 

iftxl(i.j)  .ge.  127)  then 
xl(i,j)  =  xl(i,j)-256 
else 

xl(i,j)  =  xl(i.j) 
endif 

p(i)=xl(i,j) 

80  continue 

\vritc(13'j)  p 
70  continue 

return 
end 

^  ❖  ❖  ☆  ❖  >Js  ))t  ijt  ijt  ^  ^  ^  <{s  jjt  jJe  ^  >Ji  ^  ij:  >Jc  ^  ^  j{;  ^ 

subroutine  htl(ip,jrt,sc,isx,isy,itliIiro,drow,rou'0,x0,y0) 

dimension  ip(isx.isy) 
dimension  jrt(ith.iro) 
real  sc(ith,2) 
real  x0,y0 
real  pai 

pai=  3.1415926 
jsx=isx 
jsy=isy 
jth  =  ith 
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jro  =  iro 
do  100  j=  l,jro 
do  100  i=  l,jth 
jrt(i,j)  =  0 
100  continue 

dth  =  pai/float(jth) 
do  10  i=  l,jth 
th=  float(i-l)*dth 
sc(i,l)  =  cos(th) 
sc(i,2)  =  sin(th) 

10  continue 

rowm=  sqrt(float(jsx*jsx)  +  float(jsy  *  j  sy )) 

rowO  =  rowm/2.0 

drow=  rowm/float(jro) 

x0  =  float((jsx  +  2)/2) 

y0=noat((jsy+2)/2) 

do  40  iy=  l,jsy 

y=iy 

do  30  ix=  l,jsx 
ipd  =  ip(ix,iy) 
if(ipd.le.0)  go  to  30 
x=  ix 

do  20  it=  l,jth 

row  =  (x-x0)*sc(  it,  1 )  +  (y0-y)*sc(it,2) 
r  =  (row+  row0)/drow 
ir  f=  jnint(r) 

jrt(it,ir)  =  jrt(it,ir)  +  ipd 
20  continue 

30  continue 
40  continue 
return 
end 


•>  v  £  »Ji  yf.  >;•  &  >/■  )|<  £  fy  fy  «Js  %  iiX  iji  if.  >Jc  >Ji  $  >'<  &  £  >}c  >>  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ajt  ft  ft  if.  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft 
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y.  y.  y.  y-  y 


»-  y.  y.  y.  y.  y.  y.  a  j. 


subroutine  xx 
byte  a(256) 
real  iiR(256,4096) 
integer  iiv,k,m,ic 

integer  ip(256,256),jrt(256,256),io,it,ir,ipd,ia(256) 
real  ta(256),tb(256),max,min 

real  sc(2, 256), dth,th,rowm,row0,drow,ib(256), factor, fa  1  ,fa2 

integer^  isx/256/,isy/256/,ith/256/,iro/256/ 

open(unit=  1 3,  name  =  'rcbndast.dat',  status  =  'old',  access  =  'direct' 

1  ,record$ize=  64) 

opcn(unit=  14,name=  'rcbmclu.dat', status  =  'new', access  =  'direct' 

1  ,recordsize=  64) 

c  ****read  in  the  image  data  &  convert  to  integer. *##**#**#**!**#*** 
do  10  j=  1,256 
rcod(3'j)  a 
do  20  i=  1,256 
ip(i,j)=a(i) 
if(ip(i,j).lt.0)  then 
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A 


ip(ij)=ip(ij)+256 
end  if 

20  continue 
10  continue 

maxi  =  0 
mini  =0 
do  21  j  =  1,256 
do  22  i=  1,256 

ifTip(i,j)  -gt.maxl)  maxi  =ip(i,j) 
if(ip(i,j)  .It.  mini)  mini  =  ip(i,j) 

22  continue 

21  continue 

print  900, maxi, mini 
900  format(lx,2i7) 

^  *  £  *  ffi  ^  j-  ^  5  [q  j*  iyi  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  * 

call  ht2(ip,jrt,sc,isx,isy,ith,iro,ic) 

q  tr-  #  >1"'  <s  ijc  *  <«  ^  << *  i):  j)-.  ^  #  #  ❖  <<  *  £  $  £  ijc  &  >?:  &  ft  0:  >Jt  &  #  *  * »J;  ❖  »X  fc  £  ❖  £  ❖  ❖ 

ma.x  =  jrt(  1,1) 
min  =  jrt(l,l) 
do  30  j=  1,256 
do  40  i=  1,256 

illjrt(i.j)  .gt.max)  max  =  jrt(i,j) 
if(jrt(i,j)  .It.  min)  min= jrt(i,j) 

40  continue 
30  continue 

factor = 255.0/(max-min) 
do  50  j=  1,256 
do  60  \  -  1,256 

jrt(i,j)  =  jnint((jrt(i,j)-min)*factor) 

60  continue 
50  continue 

threshold  =  255*0.2 
do  81  j=  1,256 

do  82  i=  1,256 

if(jrt(i,j).lt. threshold)  jrt(i,j)  =  0 
82  continue 

81  continue 

do  70  j=  1,256 
do  80  i=  1,256 

ifXjrt(i.j)  .ge.  127)  then 
a(i)  =  jrt(i,j)-255 
else 

a(i)  =  jrt(i.j) 


endif 

80 

continue 

\vritc(4'j)  a 

70 

continue 

•*  *  £  <t  ij.  a.  »>  *  *  *  £  £  fy  £  £  £  fc  >jt  *  *  *  £  £  *>  £  $  £  *  *  *  ❖  *!»  ❖  *  *  ❖  *  *  *  *  ❖  *  £  ❖  *  ❖  #  ❖  *  *  *  *>  *  £  £  ❖  *  $  £  ❖  *  *>  #  *  Jt  *  ❖ 

*  *  o  £  *  O  <5  *  *  £  0  £  ❖  ❖  *  *  *  £  *  »>  fc  *>  *  fc  *  <t  »!j  £  *  «jc  #  *  ijt  »ji  £  #  jjs  *  *  ifc  #  *  if  #  £  £  j£t  *  fc  £  *  i>  *  #  <t  *  *  *  *  *  *  * 


subroutine  clustering 
byte  a(256),p(256) 


71 


integer  zz(S),split/l  2S/,spIit  1  /2/,w,z,t(256) 
real  threshold 

integer  k,m,c,u(4),max,max  1  ,ta(256),tc(256),x  1  (4,4) 

integer  max2,min2,max3,isy  1/256/, ct 

integer  ip(256,256),jrt(256,256),io,it,ir,ipd,td(128) 

real  mean, mean  1  ,s(256,256),te  1 

integer  ite,tf,ib,aa(256,256),cou 

integer  adx,ady,adt,adu,ipl(4,4),irtl(4,4) 

real  tb(256),rovv(128),x(4,4) 

real  sc(4, 2), dth,th,ro\vm,ro\vO,drow, factor 

integer^  isx/2/,isy/2/,ith/2/,iro/2/ 

open(unit=  14, name = '  rcbmclu.dat', status  =  'old',  access  =  'direct' 

1  ,rccordsizc=  64) 

opcn(unit=  15,name=  'clubm.dat',status=  'new', access  =  'direct' 

1  ,recordsize  =  64) 

:  ****read  in  the  image  data  &.  convert  to  integer. ###*****##*#*###* 

cou  =  0 

do  10  j  =  l.isyl 
read(4'j)  a 
do  20  i=  l.isyl 
ip(i,j)  =  a(i) 
il(ip(i,j).lt.0)  then 

ip(i,j)-ip(i,j)+255 

endif 
jrt(i,j)  =  0 
t('0=>p(g) 

if(t(i).gt.O)  cou  =  coiH- 1 
20  continue 
10  continue 
print  902.COU 
)02  fornuu(ix,i6) 
do  $00  w=  1,10 
max  =  0 
min  =  0 

do  23  j  =  l.isyl 
do  24  i  =  l.isyl 

ilUp(i.j)-gt.max)  max-ip(i,j) 
ir(ip(i,j>.U.min)  mio  =  ip(i,j) 

M  continue 
23  continue 

print  901, max, min 

)0 1  format(  1  x/max = ',i6,'niin = ',i6) 

if(max.gt.min)  then 
ct=0 

do  100  j=  l.isyl 

do  101  i=  l.isyl 
if(ip(i,j).gc.max)  then 
ct  =  ct  *r  1 

:  tel  =  (i-l)*2*l$0/isyl 

ite  =  i 
tf=j  * 
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ib=ite**2  +  tP*2 


endif 


101  continue 
100  continue 

ite  —  ite 
ili«ib 

do  102  j=  1  ,isyl 

do  303  i=  l.isyl 

if(ip(i,j).gt.O)  then 
c  =  i**2  +  j**2 

aa(i,j)  =  (ite-i)**2  +  (tf-j)**2 
s(i,j)  =  (b  +  c-aa(i,j))/(b  +  c  +  aa(i,j)) 
if(s(i,j).lt.0.89)  then 
ip(i,j)  =  ip(i,j) 
else 

ip(i,i)=0 

endif 

endif 

ip(i.j)"ip(i»j) 

103  continue 

102  continue 
jrt(itc,tf)  =  255 
endif 


800  continue 

do  106  j  =  1  ,isy  1 

do  107  i=  l.isyl 

iRjrt(i,j).gt.I27)  then 
jrt(i,j)  =  jrt(i,j)-256 
else 

jrt(i.j)  =  jrt(i.j) 
endif 

»(0  =  jrt(i,j) 

107  continue 

\vrite(5'j)  a 
106  continue 
call  reconst 
end 

q  #  *  *  *  *  *  n«  ns  ns  s  #  ns  *  #  .is  .is  *  n<  .i:  $  «  ns  *  *  ns  *  *  *  *  n<  *  .is  is  *  *  *  *  ,i<  *  *  &  *  *  *  *  ns  *  ns  ns  ns  $  ns  ns 

£  £  >Js  >Js  £  <i  n*.  ijs  *  Is  ^  >!'•  £  >>.  *  .}&  j.v  Jji  Is  *  t*i  >Js  >Jc  &  #  ;Js  »Js  jJ;  J>.  ft  &  <<  »Js  &  *  ❖  >Ms  >};  *  ^  ❖  ❖  >Js  $  >Js  £  >Js  jfc  >Jc  £  Jjt  >Js  £  ijs  j&  >>  *  ?Jt 


subroutine  rcconst 
byte  a(256),p(256) 

integer  k,m,c,u(256),roax,maxl,ta(256),te(256),xl(256,256) 
integer  xO / 1 29/, yO / I29/,tb(256) 
integer  ip(256,256),jrt(256,256),io,it,ir,ipd,td(256) 
real  row(256),x(256,256) 

real  sc(256, 2), dth,th,rowm,rowO,drow, factor, threshold 

intcger*4  isx/256/,isy/256/,ith/256/,iro/256/ 

open(unit  =  1  ‘5, name  =  'clubm.dat', status  =  'old', access  =  'direct' 

1  ,rccordsize=  64) 

open(unit  =  1 6, name  =  'rcclubm.dat', status = 'new', access  —  'direct' 
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1  ,recordsize=  64) 

c  ****read  in  the  image  data  &  convert  to  integer. ***************** 
dro\v=2**(0.5) 
rowO=  128*drow 
do  10  j=  I  ,isy 
read(5'j)  a 
do  20  i=  l,isy 
jrt(i,j)=a(i) 

20  continue 
10  continue 

do  21  j=  l,isy 
do  22  i=  l,isy 

if(irt(i,j).lt.O)  jrt(i.j) = jrt(i,j)  +  256 
22  continue 

21  continue 

q  * *  *  * .1*.  *  *******  reconstruction  *  *  *******  *  *  *  *  *  *********  *  * 

Q  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  »Jl  **************  * 

pai  =  3.1415926 
c  =  0 

do  90  j  =  l,isy 

do  100  i=  l.isy 

if(jrt(i;j)  .gt.  0)  then 
c=c+  1 
tb(c)=  i 

tc(c)  =  j 

cndif 

100  continue 

90  continue 
print  9 1  ,c 

91  format(  lx, 'counter  =  ',i5) 
do  110  j=  l,c 

row(j)  =  tc(j)*drow-row0 
110  continue 

print  92,row(l) 

92  format(lx,'ro'vl  =  ',f6.1) 
print  93,ro\v(2) 

93  format(lx,'row2  =  ',f6.1) 
do  120  j=  l,c 

do  130  i=  l,isy 

x(i,j)  =  xO  +  (row(j)-(y0-i)*sin(tb(j)*pai/256))/cos(tb(j)’i‘pai/256) 
xl(i,j)=jnint(x(i,j)) 


td(i)=xl(i,j) 

130 

continue 

120 

continue 

do  140  j=  l,isy 

do  150  1=  l,c 

u(l)  =  xl(j,l) 

150 

continue 
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do  160  m=  l,c 
do  161  k=  l.isy 

if(k.eq.u(m))  xl(k,j)  =  255 

161  continue 

160  continue 

140  continue 

print  151,u(l),u(2) 

151  format(lx,i5,lx,i5) 
do  162  j  =  l,isy 
do  163  i=  l,isy 

if(xl(i,j).ne.255)  xl(i,j)  =  0 
163  continue 

162  continue 

do  70  j  =  l,isy 
do  80  i  =  l,isy 

il\xl(i,j)  .ge.  127)  then 
Xl(i,j)  —  xl(i,j)-256 
else 

xl(i,j)=xl(i,j) 

endif 

p(i)  =  xl(i,j) 


80 

continue 

\vrite(6'j)  p 

70 

continue 

return 

end 
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